library(rnaturalearth)
library(rnaturalearthdata)

library(sf)
library(ggplot2)
library(grid)
library(dplyr)
library(ggrepel)
library(ggtext)

set.seed(2023)


worldcoasts <- ne_coastline(scale = 110, returnclass = "sf")
countrycentroid <- read.csv("countries_centroids.csv")

### CARIBBEAN
carib_set <- data.frame(
  country = c("Antigua and Barbuda","Anguilla (UK)","Aruba (Netherlands)","Barbados",
              "Bonaire, Sint Eustatius & Saba (Netherlands)","Bahamas","Bay Islands (Honduras)",
              "Belize (Cayes)","Cayman Islands (UK) ","Dominica", "Dominican Republic","Curacao (Netherlands)",
              "Grenada","Guadeloupe (France)","Haiti","Jamaica","Martinique (France)",
              "Puerto Rico (US)","St. Kitts and Nevis","Saint Lucia",
              "Saint Martin (France)","San Andres and Providencia (Colombia)",
              "St. Vincent and the Grenadines","Trinidad and Tobago","Turks and Caicos",
              "US Virgin Islands (US)", "Virgin Islands (UK)"), 
  ISO = c("AG", "AI","AW","BB",
          "BQ","BS","HN",
          "BZ","KY","DM","DO","CW",
          "GD","GP","HT","JM","MQ",
          "PR","KN","LC",
          "MF","CO",
          "VC","TT","TC", "VI", "VG"))


### ATLANTIC
atlantic_set <- data.frame(
  country = c("Sao Tome and Principe","Cape Verde", "Bermuda (UK)"), 
  ISO = c("ST","CV","BM"))

### INDIAN
indian_set <- data.frame(
  country = c("Comoros","Maldives","Mauritius","Mayotte (France)",
              "Reunion (France)","Seychelles"), 
  ISO = c("KM","MV","MU","YT","RE","SC"))


### OCEANIA
pacific_set <- data.frame(
  country = c("American Samoa (US)","Cook Islands","Federated States of Micronesia", "Fiji",
              "French Polynesia (France)","Guam (US)","Kiribati","Marshall Islands",
              "Northern Mariana Islands (US)","New Caledonia (France)","Nauru","Papua New Guinea","Palau",
              "Samoa","Solomon Islands","Tonga", "Tuvalu", "Vanuatu"), 
  ISO = c("AS","CK","FM", "FJ",
          "PF","GU","KI","MH",
          "MP","NC", "NR","PG","PW",
          "WS","SB","TO", "TV", "VU"))


## create caribbean map

carib_set <- left_join(carib_set, countrycentroid) # merge centroid positions
carib_set <- carib_set[-c(6,7),]

#custom manual centroid positions for Belize, Honduras and Columbia (subregion being sampled, not overall country)
carib_set$longitude[which(carib_set$ISO=="BZ")] <- -88.0277
carib_set$latitude[which(carib_set$ISO=="BZ")] <- 17.7612

carib_set$longitude[which(carib_set$ISO=="HN")] <- -85.8793
carib_set$latitude[which(carib_set$ISO=="HN")] <- 16.4827

carib_set$longitude[which(carib_set$ISO=="CO")] <- -81.7051
carib_set$latitude[which(carib_set$ISO=="CO")] <- 12.5769

set.seed(2003)
#  coastline basemap
carib <- ggplot() +
  geom_sf(data = worldcoasts) +
  coord_sf(xlim = c(-90, -50), ylim = c(9,29)) +
  theme_bw()+
  theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank(), #remove x axis ticks
        axis.text.y=element_blank(),  #remove y axis labels
        axis.ticks.y=element_blank(),  #remove y axis ticks
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()
  )

#  generate ggrepel positions w/ manual adjustments
set.seed(2023)

carib_xnudge <- c(10, #AG
                  7, #AI
                  5, #AW
                  3, #BB
                  5, #BQ
                  5, #BS
                  -2, #HN
                  -10, #BZ
                  0, #KY
                  8, #DM
                  0, #DO
                  0, #CW
                  5, #GD
                  6, #GP
                  -4,#HT
                  -5,#JM
                  8,#MQ
                  -2, #PR
                  9,#KN
                  9,#LC
                  8,#MF
                  0,#CO
                  7,#VC
                  8,#TT
                  3,#TC
                  6,#VI
                  4)#VG

carib_ynudge <- c(10, #AG
                  8, #AI
                  -2, #AW
                  -6, #BB
                  -3, #BQ
                  0, #BS
                  0, #HN
                  0, #BZ
                  0, #KY
                  0, #DM
                  2, #DO
                  -1, #CW
                  -5, #GD
                  -2, #GP
                  5,#HT
                  -5,#JM
                  0,#MQ
                  -2, #PR
                  2,#KN
                  -2,#LC
                  0,#MF
                  0,#CO
                  0,#VC
                  -2,#TT
                  5,#TC
                  5,#VI
                  10)#VG


p1 <- carib + geom_text_repel(data = carib_set, aes(x = longitude, y = latitude,label=ISO),
                              box.padding = unit(2, "lines"), size=6,
                              position = position_nudge_repel(x=carib_xnudge, y=carib_ynudge)) +
  geom_point(data = carib_set, aes(x = longitude, y = latitude))

quartz(width=7, height=5.5)

# Get x and y plot ranges 
set.seed(2023)
ggp1 <- ggplot_build(p1)
xrg <- c(-90,-50)
yrg <- c(9,29)
set.seed(2023)
p1
set.seed(2023)
grid.force()
kids <- childNames(grid.get("textrepeltree", grep = TRUE))
kids <- kids[grep("textrepelgrob",kids)]


# Function: get the x and y positions of a single ggrepel label
get.xy.pos.labs <- function(n) {
  set.seed(2023)
  grb <- grid.get(n)
  data.frame(
    x = xrg[1]+diff(xrg)*convertX(grb$x, "native", valueOnly = TRUE),
    y = yrg[1]+diff(yrg)*convertY(grb$y, "native", valueOnly = TRUE)
  )
}

# Get positions of all ggrepel labels
set.seed(2023)
dts <- do.call(rbind, lapply(kids, get.xy.pos.labs))
dts$ISO <- carib_set$ISO

carib_set2 <- left_join(carib_set, dts)

carib_set2$label <- paste0("<img src='/map_figures/",carib_set2$ISO,"_map.png' width='50'/>") # Main Text Figure 2
#carib_set2$label <- paste0("<img src='/map_figures/",carib_set2$ISO,"_map_climate.png' width='50'/>") # SI Figure 3

### FINAL MAP

carib_set2$x_mod <- carib_set2$x
carib_set2$y_mod <- carib_set2$y

carib_set2$y_mod[which(carib_set2$ISO=="BS")] <- carib_set2$y[which(carib_set2$ISO=="BS")]+0.75

carib_set2$x_mod[which(carib_set2$ISO=="TC")] <- carib_set2$x[which(carib_set2$ISO=="TC")]-4
carib_set2$y_mod[which(carib_set2$ISO=="TC")] <- carib_set2$y[which(carib_set2$ISO=="TC")]-2


carib_set2$x_mod[which(carib_set2$ISO=="PR")] <- carib_set2$x[which(carib_set2$ISO=="PR")]+7

carib_set2$x_mod[which(carib_set2$ISO=="VG")] <- carib_set2$x[which(carib_set2$ISO=="VG")]+3
carib_set2$y_mod[which(carib_set2$ISO=="VG")] <- carib_set2$y[which(carib_set2$ISO=="VG")]+5

carib_set2$x_mod[which(carib_set2$ISO=="VI")] <- carib_set2$x[which(carib_set2$ISO=="VI")]-5
carib_set2$y_mod[which(carib_set2$ISO=="VI")] <- carib_set2$y[which(carib_set2$ISO=="VI")]+3

carib_set2$y_mod[which(carib_set2$ISO=="MF")] <- carib_set2$y[which(carib_set2$ISO=="MF")]+4

carib_set2$x_mod[which(carib_set2$ISO=="AG")] <- carib_set2$x[which(carib_set2$ISO=="AG")]+7
carib_set2$y_mod[which(carib_set2$ISO=="AG")] <- carib_set2$y[which(carib_set2$ISO=="AG")]+7

carib_set2$x_mod[which(carib_set2$ISO=="KN")] <- carib_set2$x[which(carib_set2$ISO=="KN")] -4.5
carib_set2$y_mod[which(carib_set2$ISO=="KN")] <- carib_set2$y[which(carib_set2$ISO=="KN")] -1

carib_set2$x_mod[which(carib_set2$ISO=="DO")] <- carib_set2$x[which(carib_set2$ISO=="DO")]
carib_set2$y_mod[which(carib_set2$ISO=="DO")] <- carib_set2$y[which(carib_set2$ISO=="DO")]+2

carib_set2$x_mod[which(carib_set2$ISO=="JM")] <- carib_set2$x[which(carib_set2$ISO=="JM")]+3
carib_set2$y_mod[which(carib_set2$ISO=="JM")] <- carib_set2$y[which(carib_set2$ISO=="JM")]+2

carib_set2$y_mod[which(carib_set2$ISO=="HN")] <- carib_set2$y[which(carib_set2$ISO=="HN")]-2

carib_set2$y_mod[which(carib_set2$ISO=="CO")] <- carib_set2$y[which(carib_set2$ISO=="CO")]-2
carib_set2$x_mod[which(carib_set2$ISO=="CO")] <- carib_set2$x[which(carib_set2$ISO=="CO")]-2

carib_set2$x_mod[which(carib_set2$ISO=="KY")] <- carib_set2$x[which(carib_set2$ISO=="KY")]-3
carib_set2$y_mod[which(carib_set2$ISO=="KY")] <- carib_set2$y[which(carib_set2$ISO=="KY")]+3

carib_set2$x_mod[which(carib_set2$ISO=="BZ")] <- carib_set2$x[which(carib_set2$ISO=="BZ")]+3
carib_set2$y_mod[which(carib_set2$ISO=="BZ")] <- carib_set2$y[which(carib_set2$ISO=="BZ")]-6

carib_set2$x_mod[which(carib_set2$ISO=="AW")] <- carib_set2$x[which(carib_set2$ISO=="AW")]-1
carib_set2$y_mod[which(carib_set2$ISO=="AW")] <- carib_set2$y[which(carib_set2$ISO=="AW")]+3

carib_set2$y_mod[which(carib_set2$ISO=="DM")] <- carib_set2$y[which(carib_set2$ISO=="DM")]+1

carib_set2$x_mod[which(carib_set2$ISO=="MQ")] <- carib_set2$x[which(carib_set2$ISO=="MQ")]+5
carib_set2$y_mod[which(carib_set2$ISO=="MQ")] <- carib_set2$y[which(carib_set2$ISO=="MQ")]+2.5

carib_set2$x_mod[which(carib_set2$ISO=="CW")] <- carib_set2$x[which(carib_set2$ISO=="CW")]-3
carib_set2$y_mod[which(carib_set2$ISO=="CW")] <- carib_set2$y[which(carib_set2$ISO=="CW")]-2

carib_set2$y_mod[which(carib_set2$ISO=="BQ")] <- carib_set2$y[which(carib_set2$ISO=="BQ")]-3

carib_set2$y_mod[which(carib_set2$ISO=="TT")] <- carib_set2$y[which(carib_set2$ISO=="TT")]-2

carib_set2$y_mod[which(carib_set2$ISO=="BB")] <- carib_set2$y[which(carib_set2$ISO=="BB")]-1

carib_set2$x_mod[which(carib_set2$ISO=="VC")] <- carib_set2$x[which(carib_set2$ISO=="VC")]+3
carib_set2$y_mod[which(carib_set2$ISO=="VC")] <- carib_set2$y[which(carib_set2$ISO=="VC")]+4

carib_set2$x_mod[which(carib_set2$ISO=="GD")] <- carib_set2$x[which(carib_set2$ISO=="GD")]+3

set.seed(2023)

carib +
  geom_point(data = carib_set, aes(x = longitude, y = latitude))+
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank())+
  geom_segment(data=carib_set2, aes(x=longitude, y=latitude, xend=x_mod, yend=y_mod))+
  geom_richtext(data=carib_set2, aes(x=x_mod, y=y_mod, label=label),
                label.padding=grid::unit(rep(0, 4), "pt"),
                label.r=unit(0, "lines"),
                fill = NA, label.color = NA)



## create Atlantic map

atlantic_set <- left_join(atlantic_set, countrycentroid) # merge centroid positions


#  coastline basemap
atlantic <- ggplot() +
  geom_sf(data = worldcoasts) +
  coord_sf(xlim = c(-65, 5), ylim = c(0,35)) +
  theme_bw()+
  theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank(), #remove x axis ticks
        axis.text.y=element_blank(),  #remove y axis labels
        axis.ticks.y=element_blank(),  #remove y axis ticks
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()
  )

#  generate ggrepel positions w/ manual adjustments
set.seed(2023)

p2 <- atlantic + geom_text_repel(data = atlantic_set, aes(x = longitude, y = latitude,label=ISO),
                                 box.padding = unit(2, "lines"), size=6) +
  geom_point(data = atlantic_set, aes(x = longitude, y = latitude))


# Get x and y plot ranges 
set.seed(2023)
ggp2 <- ggplot_build(p2)
xrg <- c(-65, 5)
yrg <- c(0,35)
set.seed(2023)
p2
set.seed(2023)
grid.force()
kids <- childNames(grid.get("textrepeltree", grep = TRUE))
kids <- kids[grep("textrepelgrob",kids)]


# Function: get the x and y positions of a single ggrepel label
get.xy.pos.labs <- function(n) {
  set.seed(2023)
  grb <- grid.get(n)
  data.frame(
    x = xrg[1]+diff(xrg)*convertX(grb$x, "native", valueOnly = TRUE),
    y = yrg[1]+diff(yrg)*convertY(grb$y, "native", valueOnly = TRUE)
  )
}

# Get positions of all ggrepel labels
set.seed(2023)
dts <- do.call(rbind, lapply(kids, get.xy.pos.labs))
dts$ISO <- atlantic_set$ISO

atlantic_set2 <- left_join(atlantic_set, dts)

atlantic_set2$label <- paste0("<img src='/map_figures/",atlantic_set2$ISO,"_map.png' width='50'/>") # Main Text Figure 2
#atlantic_set2$label <- paste0("<img src='/map_figures/",atlantic_set2$ISO,"_map_climate.png' width='50'/>") # SI Figure 3

### FINAL MAP

atlantic_set2$x_mod <- atlantic_set2$x
atlantic_set2$y_mod <- atlantic_set2$y

set.seed(2023)

quartz(width=3.5, height=2)
atlantic +
  geom_point(data = atlantic_set, aes(x = longitude, y = latitude))+
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank())+
  geom_segment(data=atlantic_set2, aes(x=longitude, y=latitude, xend=x_mod, yend=y_mod))+
  geom_richtext(data=atlantic_set2, aes(x=x_mod, y=y_mod, label=label),
                label.padding=grid::unit(rep(0, 4), "pt"),
                label.r=unit(0, "lines"),
                fill = NA, label.color = NA)



# create Indian Map
indian_set <- data.frame(
  country = c("Comoros","Maldives","Mauritius","Mayotte (France)",
              "Reunion (France)","Seychelles"), 
  ISO = c("KM","MV","MU","YT","RE","SC"))

indian_set <- left_join(indian_set, countrycentroid) # merge centroid positions


#  coastline basemap
indian <- ggplot() +
  geom_sf(data = worldcoasts) +
  coord_sf(ylim = c(-24, 1), xlim = c(40,75)) +
  theme_bw()+
  theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank(), #remove x axis ticks
        axis.text.y=element_blank(),  #remove y axis labels
        axis.ticks.y=element_blank(),  #remove y axis ticks
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()
  )

#  generate ggrepel positions w/ manual adjustments
set.seed(2023)

p3 <- indian + geom_text_repel(data = indian_set, aes(x = longitude, y = latitude,label=ISO),
                               box.padding = unit(2, "lines"), size=6) +
  geom_point(data = indian_set, aes(x = longitude, y = latitude))


# Get x and y plot ranges 
set.seed(2023)
ggp3 <- ggplot_build(p3)
xrg <- c(40,75)
yrg <- c(-24,1)
set.seed(2023)
p3
set.seed(2023)
grid.force()
kids <- childNames(grid.get("textrepeltree", grep = TRUE))
kids <- kids[grep("textrepelgrob",kids)]


# Function: get the x and y positions of a single ggrepel label
get.xy.pos.labs <- function(n) {
  set.seed(2023)
  grb <- grid.get(n)
  data.frame(
    x = xrg[1]+diff(xrg)*convertX(grb$x, "native", valueOnly = TRUE),
    y = yrg[1]+diff(yrg)*convertY(grb$y, "native", valueOnly = TRUE)
  )
}

# Get positions of all ggrepel labels
set.seed(2023)
dts <- do.call(rbind, lapply(kids, get.xy.pos.labs))
dts$ISO <- indian_set$ISO

indian_set2 <- left_join(indian_set, dts)

indian_set2$label <- paste0("<img src='/map_figures/",indian_set2$ISO,"_map.png' width='50'/>") # Main Text Figure 2
#indian_set2$label <- paste0("<img src='/map_figures/",indian_set2$ISO,"_map_climate.png' width='50'/>") ## SI Figure 3

### FINAL MAP

indian_set2$x_mod <- indian_set2$x
indian_set2$y_mod <- indian_set2$y

indian_set2$y_mod[which(indian_set2$ISO=="RE")] <- indian_set2$y[which(indian_set2$ISO=="RE")]-6
indian_set2$x_mod[which(indian_set2$ISO=="RE")] <- indian_set2$x[which(indian_set2$ISO=="RE")]+13

indian_set2$y_mod[which(indian_set2$ISO=="YT")] <- indian_set2$y[which(indian_set2$ISO=="YT")]-2
indian_set2$x_mod[which(indian_set2$ISO=="YT")] <- indian_set2$x[which(indian_set2$ISO=="YT")]+2

set.seed(2023)



quartz(width=3.5, height=2)
indian +
  geom_point(data = indian_set, aes(x = longitude, y = latitude))+
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank())+
  geom_segment(data=indian_set2, aes(x=longitude, y=latitude, xend=x_mod, yend=y_mod))+
  geom_richtext(data=indian_set2, aes(x=x_mod, y=y_mod, label=label),
                label.padding=grid::unit(rep(0, 4), "pt"),
                label.r=unit(0, "lines"),
                fill = NA, label.color = NA)




### OCEANIA
pacific_set <- data.frame(
  country = c("American Samoa (US)","Cook Islands","Federated States of Micronesia", "Fiji",
              "French Polynesia (France)","Guam (US)","Kiribati","Marshall Islands",
              "Northern Mariana Islands (US)","New Caledonia (France)","Nauru","Papua New Guinea","Palau",
              "Samoa","Solomon Islands","Tonga", "Tuvalu", "Vanuatu"), 
  ISO = c("AS","CK","FM", "FJ",
          "PF","GU","KI","MH",
          "MP","NC", "NR","PG","PW",
          "WS","SB","TO", "TV", "VU"))



# create Pacific Map
quartz(width=7, height=5.5)

pacific_set <- left_join(pacific_set, countrycentroid) # merge centroid positions
pacific_set$longitude_mod <- pacific_set$longitude
pacific_set$longitude_mod[which(pacific_set$longitude<0)] <- pacific_set$longitude_mod[which(pacific_set$longitude<0)]+360

#  coastline basemap
pacific <- ggplot() +
  geom_sf(data = worldcoasts) +
  coord_sf(xlim = c(130, 220), ylim = c(-25,20)) +
  theme_bw()+
  theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank(), #remove x axis ticks
        axis.text.y=element_blank(),  #remove y axis labels
        axis.ticks.y=element_blank(),  #remove y axis ticks
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()
  )

#  generate ggrepel positions w/ manual adjustments
set.seed(2023)

p4 <- pacific + geom_text_repel(data = pacific_set, aes(x = longitude_mod, y = latitude,label=ISO),
                                box.padding = unit(2, "lines"), size=6) +
  geom_point(data = pacific_set, aes(x = longitude_mod, y = latitude))


# Get x and y plot ranges 
set.seed(2023)
ggp4 <- ggplot_build(p4)
xrg <- c(130,220)
yrg <- c(-25,20)
set.seed(2023)
p4
set.seed(2023)
grid.force()
kids <- childNames(grid.get("textrepeltree", grep = TRUE))
kids <- kids[grep("textrepelgrob",kids)]


# Function: get the x and y positions of a single ggrepel label
get.xy.pos.labs <- function(n) {
  set.seed(2023)
  grb <- grid.get(n)
  data.frame(
    x = xrg[1]+diff(xrg)*convertX(grb$x, "native", valueOnly = TRUE),
    y = yrg[1]+diff(yrg)*convertY(grb$y, "native", valueOnly = TRUE)
  )
}

# Get positions of all ggrepel labels
set.seed(2023)
dts <- do.call(rbind, lapply(kids, get.xy.pos.labs))
dts$ISO <- pacific_set$ISO

pacific_set2 <- left_join(pacific_set, dts)

pacific_set2$label <- paste0("<img src='/map_figures/",pacific_set2$ISO,"_map.png' width='50'/>") # Main Text Figure 2
#pacific_set2$label <- paste0("<img src='/map_figures/",pacific_set2$ISO,"_map_climate.png' width='50'/>") # SI Figure 3 


### FINAL MAP

pacific_set2$x_mod <- pacific_set2$x
pacific_set2$y_mod <- pacific_set2$y

pacific_set2$y_mod[which(pacific_set2$ISO=="PW")] <- pacific_set2$y[which(pacific_set2$ISO=="PW")]+5
pacific_set2$x_mod[which(pacific_set2$ISO=="PW")] <- pacific_set2$x[which(pacific_set2$ISO=="PW")]

pacific_set2$y_mod[which(pacific_set2$ISO=="PG")] <- pacific_set2$y[which(pacific_set2$ISO=="PG")]+5
pacific_set2$x_mod[which(pacific_set2$ISO=="PG")] <- pacific_set2$x[which(pacific_set2$ISO=="PG")]

set.seed(2023)

pacific_set <- pacific_set[-which(pacific_set$ISO=="CK"),] # no survey data
pacific_set2 <- pacific_set2[-which(pacific_set2$ISO=="CK"),] # no survey data

quartz(width=7, height=5.5)
pacific +
  geom_point(data = pacific_set, aes(x = longitude_mod, y = latitude))+
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank())+
  geom_segment(data=pacific_set2, aes(x=longitude_mod, y=latitude, xend=x_mod, yend=y_mod))+
  geom_richtext(data=pacific_set2, aes(x=x_mod, y=y_mod, label=label),
                label.padding=grid::unit(rep(0, 4), "pt"),
                label.r=unit(0, "lines"),
                fill = NA, label.color = NA)
